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ABSTRACT 

We report theoretical evidence that bulk nonlinear materials weakly interacting with highly localized plasmonic modes in 
ultra-sub-wavelength metallic nanostructures can lead to nonlinear effects at the single plasmon level in the visible range. In 
particular, the two-plasmon interaction energy in such systems is numerically estimated to be comparable with the typical 
plasmon linewidths. Localized surface plasmons are thus predicted to exhibit a purely nonclassical behavior, which can be 
clearly identified by a sub-Poissonian second-order correlation in the signal scattered from the quantized plasmonic field 
under coherent electromagnetic excitation. We explicitly show that systems sensitive to single-plasmon scattering can be 
experimentally realized by combining electromagnetic confinement in the interstitial region of gold nanodimers with local 
infiltration or deposition of ordinary nonlinear materials. We also propose configurations that could allow to realistically detect 
such an effect with state-of-the-art technology, overcoming the limitations imposed by the short plasmonic lifetime. 


Introduction 

Strongly confined electromagnetic fields at the nanoscale allow to enhance radiation-matter interaction in the optical or 
near-infrared domain to such a level that scattering of single radiation quanta may become important. On the one hand, this is a 
key target of modern nanophotonics research, e.g., for the realization of all-optical sensing, switching, and routing at the lowest 
possible control power, up to the level of a single energy quantum 1 . On the other hand, it constitutes a potential playground to 
explore fundamental quantum manybody physics related to strongly correlated states in open systems 2 . Within this context, 
single-photon nonlinear optics has been predicted in a variety of nanophotonic systems 3-10 , and demonstrated through the 
single-photon blockade effect, i.e., inhibition of the absorption of a single quantum of energy due to the presence of another 
one, both in atomic 11 and solid-state 12-14 cavity-quantum-electrodynamics (CQED). In these realizations, a single quantum 
emitter — an atom or a quantum dot — strongly interacts with a single confined mode of a suitably engineered electromagnetic 
resonator. 

Among the different routes to the electromagnetic confinement, plasmonic systems allow for enhancement of radiation- 
matter interaction to scales well below the diffraction limit, hardly reachable with conventional optical means. In particular, 
surface plasmons, i.e., mixed radiation-matter excitations that arise at metal-dielectric interfaces due to the coupled oscillation 
of the electromagnetic field and the electron charge density 15 , are currently attracting considerable interest to develop quantum 
plasmonic devices 16,17 . One of the ultimate goals is to explore the use of surface plasmon excitations to enhance nonlinear 
interactions with single quantum emitters 18-20 . From the point of view of classical electromagnetism, optical bistability in 
metal nanoantennas involving a third-order nonlinear medium has been extensively analyzed 21 . On the other hand, in the limit 
of extreme localization of the plasmonic field, a truly quantum regime can be reached, characterized by the nonlinear interaction 
among single quanta of energy. Only a few pioneering works have proposed to reach such regime of single-plasmon blockade 
by using non-resonant material nonlinearities. In particular, early attempts of measuring a plasmonic analog of the Coulomb 
blockade effect have been reported in Ref . 22 by detecting quantized steps in the low-power transmission through sub-wavelength 
metallic pinholes. More recently, theoretical proposals were focusing on confined surface plasmons in monolayer graphene 
nanocavities 2 - 1,24 , to realize a plasmon blockade effect in the mid-infrared range. As continuous-wave excitation is generally 
assumed, the actual observability of such an effect in a realistic experimental set-up may be strongly hindered by the extremely 
short plasmonic lifetimes, especially when turning to conventional noble metals at optical frequencies. In general, there has 
been little consideration on the possibility of overcoming these limitations with pulsed excitation techniques. 

At variance from the conventional CQED systems, which typically rely on long lifetimes, thus limiting the maximum 
single-photon emission rates, here we propose a system where the intrinsically small material nonlinearity 25 can be enhanced 



Figure 1 . (a) Sketch of the system under consideration in this work. A nanodimer made up of two metallic nanoparticles 
(MNPs) is embedded in a nonlinear medium with a third-order susceptibility and it is coherently pumped at frequency (0 p 
from the exterior. The dimer sustains a localized surface plasmon with a decay rate yo = y ra( j + y nr , due to radiation losses and 
dissipation, (b) Scheme of the energy levels of the system, highlighting the single-plasmon-blockade mechanism. 


due to electromagnetic field confinement, so that a regime of quantum nonlinearity can be reached even without the use 
of long-lifetime quantum emitters. This potentially provides an extremely fast source of quantum states of radiation with 
order-of-magnitude improvements over full-dielectric systems, to be eventually employed in quantum information. Specifically, 
in the present work we analyze the quantized surface plasmonic excitations of a metal nanodimer, a nanostructure largely 
within reach of state-of-the-art technology in terms of size and gaps between its constitutive elements 26-29 . Localized surface 
plasmons are characterized by a finite linewidth due to both radiation losses and dissipation inside the metal. We show that 
the single-plasmon blockade can be reached quite straightforwardly under an external coherent excitation, e.g., an external 
laser input, by infiltrating the interstitial region of the dimer with a sufficiently nonlinear optical material (e.g., molecular dye 
solutions), as sketched in Fig. 1(a). Even if a similar system was proposed in Ref. 22 , the actual probing of plasmon blockade in 
terms of quantum correlation measurements, which we theoretically address in the present work, was not considered before. 

Although it is easy to compute the response of the metallic nanoparticle with different kinds of external excitations, such as 
plane waves or oscillating dipoles, establishing a quantum-mechanical theory of the plasmonic excitations is a far more complex 
task, due to the intrinsically lossy (radiative and non-radiative) character of surface plasmons 21 ’■ 24 ' 40 - 42 . In the following, we 
show that, by relying on the formalism of quasinormal modes, it is indeed possible to rigorously demonstrate that the system 
can be modeled with a quantum master equation involving a single bosonic operator, similarly to non-absorbing photonic 
systems. Within this theoretical picture, the physical basis of the effect can be grasped by looking at the first excitation levels of 
the system, which are schematically reported in Fig. 1(b). We suppose to coherently pump the nanodimer at a frequency co p , 
in resonance with the bare surface plasmon energy Picoq. The simultaneous presence of two highly-confined plasmons in the 
spatial region occupied by the nonlinear material is associated with a large nonlinear interaction, U, which produces a shift of 
the double-excitation energy. The magnitude of U essentially depends on the degree of localization of the electromagnetic 
field, which can be quantified through the effective volume of the plasmonic mode, V e ff. As we discuss in the following, 
although the formal expression for 14 rt is the same as for dielectric systems, it is essential to ensure the correct normalization 
of the plasmon eigenfield in agreement with the theory of quasinormal modes. Moreover, due to their intrinsically lossy 
character, plasmonic levels present a finite linewidth, yo- When the effective volume is sufficiently small (for the system 
under consideration, jS 10 2 /im 3 ), the magnitude of U could result in a shift of the two-plasmon excitation frequency 
(A co = 2U / h , see the scheme in Fig. 1) that is larger than the natural linewidth of the mode, calculated to be around 100 meV 
for the typical nanodimers under consideration. As a result, the system cannot absorb a second plasmon but upon re-emission 
of the first one, thus effectively becoming a source of single plasmonic excitations. 

Finally, we theoretically show that single-plasmon blockade could be experimentally measured, despite the fs-scale plasmon 
lifetime, by using a pulsed excitation source. Besides opening another route to the experimental study of quantum plasmonic 
effects, truly meant as the mutual interaction of single plasmon excitations at the nanoscale, these results could enable an 
ultra-high-rate source of single radiation quanta at visible wavelengths. 
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Figure 2. (a) Normalized decay rate of a dipolar emitter located at the center of the 10-nm interstitial region between two 
60-nm-diameter, 15-nm-high gold nanodisks (see inset), v,v the emitter energy. The calculation is performed with the BEM, 
DDA, and Fermi golden rule (FGR) applied to Eq. (4). (b) Cross section of the intensity of the QNM electric field, A| 2 , along 
the dimer middle plane. The units for the QNM field are chosen in agreement with Eq. (4). 


Results and discussion 

We consider localized surface plasmons of gold nanodimers in the context of the boundary element method (BEM) formalism 
developed in Ref. 30, which can be applied to any system of locally homogeneous regions of space separated by abrupt interfaces. 
In Fig. 2(a) we plot the semiclassical decay rate of a dipole emitter with momentum p located at the center of the 10-nm gap 
between two gold nanodisks and directed along the inter-particle axis (see inset), showing a clear peak around ft) = 1.7 eV due 
to coupling with a (longitudinal) surface plasmon of the nanoantenna. In this work, we use the experimental dielectric function 
of evaporated gold 31 and we assume a background refractive index = 1.5, to account in an averaged way for the surrounding 
dielectrics. The rate computed with the MNPBEM toolbox (blue solid curve), a publicly available implementation of the 
BEM 32 , is compared with a completely independent calculation (dots) within the discrete dipole approximation (DDA) 33 34 , 
showing very good agreement. Both quantities are normalized to the dipole free-space decay rate Ff ree = ft) 3 /; is/) 2 /(37T£o/ic 3 ). 

It is recognized that the definition of localized plasmons in metal nanostructures, and, more generally, of the natural 
oscillation modes of leaky optical systems, can be made rigorous in the framework of quasinormal modes (QNMs) 35-37 , i.e., 
the solutions of a non-Hermitian differential equation with complex eigenfrequencies ft). We have computed QNMs with the 
following procedure. In the BEM formalism, the electric field in each region is decomposed into an incident and a scattered 
part, E(r, o) = E; nc (r, ft)) +E sc (r,ft>), the latter being identified with the field generated by a (fictitious) charge and current 
distribution on the enclosing surface. After discretizing the surface into a collection of N representatives points, the charge and 
current distribution is calculated from the solution of a linear problem of the form Ex = a, where E(ft)) is a N xN matrix and a 
is a vector constructed from the incident field. For brevity, we omit the full equations of the method and we refer to Ref. 30 for 
the definition of the involved quantities (see Supplementary Information). For our purposes, it suffices to notice that the QNMs 
of the system are obtained by the condition detE(a>) = 0, which corresponds to the nonlinear eigenvalue problem 

E(cS)x = 0 (1) 

for the generalized eigenvector x. More importantly, the quasinormal field (QNF), A(r), can be computed as the electric field 
generated by the surface charge and current distribution obtained from the eigenvector (see Supplementary Information). We 
have solved Eq. (1) working within the MNPBEM toolbox and using the two-sided Rayleigh functional iterative algorithm' 8 . 
For instance, for the plasmon associated to the peak in Fig. 2(a), we obtain ft) ~ 1.68 — /0.07 e V, and the QNF intensity plotted 
in Fig. 2(b). 

At variance from the well-known relation f dre(r)|<?(r)| 2 = 1, which holds for normal modes in non-dispersive media, 
QNMs satisfy a much more complex normalization condition'^ - ' 7 . Instead of explicitly calculating the normalization integral, 
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we follow an implicit approach to normalize the field, similar to that proposed in Ref. 39. Indicating with ro the position of the 
dipole emitter at the center of the dimer gap, we normalize the QNM by imposing the following relation ■ for ft) —> ft): 

E sc (r, ft)) ~ -a> * /f^ r °L ^(r) ( 2 ) 

2£o (CO — co) 

(the factor 2eo is for dimensional convenience). We further improved the procedure by factoring out the singular term and 
imposing the relation only to the residues, with great advantages in terms of numerical accuracy and computational efficiency 
(see Supplementary Information). 

The usual starting point for the quantization of plasmonic systems is the system-bath approach 40 , where the electric field is 
expanded on a continuum of frequency-dependent bosonic operators. Following Ref. 40 and introducing the collective modes 
/3 (g)), Eq. (2) implies that for ft) ss 0 )q and Yur/c <C 1 the field operator can be approximated as 




d "7§ /i( " K 


H.c., 


(3) 


with g 2 (co) = 7o/[(ft) — G)o) 2 + ')o/4] and <^(r) the QNF introduced previously (see Supplementary Information). This model 
corresponds to a structured bosonic bath with a Lorentzian spectral density, whose central frequency coo = Re(fti) and linewidth 
7 o = — 2Im(cd) are directly related to the QNM eigenfrequency. According to a well-known theoretical result 40,41 , the dynamics 
of such a model is equivalent to that of a single bosonic mode, p, which can thus be interpreted as the plasmon destruction 
operator, coupled to a flat reservoir with dissipation rate yo- Thus, the electric field operator can be equivalently written 


E(r,0 =i \ 5 ’ ( 4 ) 

with the system dynamics being described by a density-matrix master equation in the Markov approximation: 

ihp = [H, p] + ^ [2/5p/3 + - /3 f /3p - p/3 T /3]. (5) 

This represents a commonly employed approach in the recent field of quantum plasmonics 20,24,40 ’ 42 . 

To check its validity, in Fig. 2 we also show (red dashed curve) the emitter decay rate calculated in the weak coupling 
regime by applying the Fermi golden rule to Eq. (4), i.e., r(o) = ft)o7o|p ■ <#(ro)| 2 /[2eo7f((ft) — coo) 2 + 7j/ 4)] and normalized 
to rfree(ft)). The very good agreement with the semi-classical result near the resonance frequency confirms the validity of our 
model and justifies the single-mode approximation of Eq. (4). 

We assume the plasmonic dimer to be embedded in a nonlinear medium, e.g., with a strong third-order nonlinear suscepti¬ 
bility %' 31 , and we write the Hamiltonian in the form 

Ho = Sft)o/3 t p-)-t//3'/3 t /3/3 (6) 

where the first term accounts for the energy of the plasmon mode, including the contribution of metal dispersion (as implied by 
the QNM normalization condition). The second term is a four-operator product deriving from the third-order susceptibility 
of the nonlinear medium around the dimer and accounting for two-plasmon scattering. In a sense, this term is the analog of 
Hubbard interaction in electron systems. If we assume the X ' 31 to be nonzero only in a region outside the nanodimers and 
dispersionless, we can use the same expression for the nonlinear interaction energy already introduced for dielectric systems 8 : 

u= (heoot [ drx (3) {r) \g ir) \\ (7) 

£o J 

The master equation approach would also allow to include nonlinear loss terms, such as the nonlinear absorption introduced 
by an imaginary part of the X' 2 ' 1 response, which we neglect due to the small pumping rates considered here (see, e.g.. Ref. 8 ). 
Moreover, we are considering a scalar nonlinear response (i.e., no tensorial coupling between the different S' components), 
since we are mainly interested in a proof-of-principle theoretical demonstration of the quantum nonlinear behavior in a standard 
experimental setup as a function of the possible values of the X ^ ' elements. In this respect, we can simplify the nonlinear 
coupling energy to U = (/ift)o) 2 %* 3 *V r eff 1 /eo, where the effective volume of the confined plasmon is defined as 

V*=J a drK(r)| 4 , (8) 
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Figure 3. (a) Effective length, V ef £ , and (b) decay rate, Tiyo, for the lowest-frequency longitudinal plasmon in dimers of 
different gold nanoparticles (see insets), versus the interstitial gap size. The curves refer to the case of disks with diameter of 
40, 60, and 80 nm and bowtie dimers obtained from equilater triangles with edge of 40 and 60 nm (curvature radius = 5 nm). 
The height of all particles is 15 nm. 


the integration being limited to the domain il represented by the bounded physical volume filled with the nonlinear medium. 
Veit essentially reflects the degree of localization of the QNF in region Q. 

As a final approximation, we are neglecting any nonlocal effects of the x^ response, which are known to arise for gap sizes 
certainly smaller than 1 nm 4 ’' 44 . We also stress that recent experimental investigations at the single or few molecules level 
within the sub-nm gap size between a metallic microsphere and a planar surface have been thoroughly explained by purely 
local theoretical modeling 40 . 

By solving the eigenproblem (1) within the MNPBEM toolbox, we have characterized the lowest-frequency longitudinal 
plasmon modes of two kinds of nanodimers: nanodisk dimers and bowtie dimers made of two mirroring triangular nanoparticles, 
for various sizes and inter-particle distances. The results for the effective length E eff and the linewidth hyo are summarized in 
Fig. 3(a) and (b), respectively. The effective length depends mostly on the gap size, reflecting the plasmonic field enhancement 
in the interstitial region, whereas the linewidth is mainly affected by the particle size (due to radiative decay) and plasmon 
frequency (due to metal dissipation). Single plasmon blockade is based on a delicate interplay between V e ff and Jq, thus both 
quantities must be taken into account when choosing the best geometry. It is already evident, however, that bowtie dimers are 
favored with respect to nanodisks. In any case, we notice that even the smallest gap size considered in our simulations, i.e. 5 
nm, is much larger than what is currently achieved 28,29 . The integral in Eq. (8) has been performed in the region external to the 
metal and bounded by a box spanning the same height of the dimer along z and extending 10 nm beyond the particle extremities 
in the xy plane. This choice is compatible, for instance, with a film of nonlinear medium coated on top of a substrate. In general, 
we observed little variation of the values of the effective length when varying the bounding box, as long as the interstitial 
region between the metal nanostructures is entirely included in the integration volume. We also notice that the structures under 
consideration are well described within the regime of classical electromagnetic theory, and do not require to include quantum 
effects due to the microscopic nature of the metallic surfaces 46,47 . 

Following our previous considerations, the key quantity for quantifying the occurrence of blockade effects in the system is 
the ratio between the double-excitation energy shift and the plasmon linewidth, or, equivalently, the ratio f//(/iyo). In Fig. 4 we 
plot this ratio as a function of the nonlinear susceptibility of the infiltrated medium for a few cases of potential interest, selected 
from the results presented in Fig. 3. Three possible nanoplasmonic dimers are considered, either with disk or bowtie geometry. 
We expect that plasmon blockade starts to play a significant role in the quantum dynamics of the system when the energy shift 
becomes comparable with the linewidth, i.e., when U/(hyo) > As it can be seen from Fig. 4, this condition corresponds to 
values of the nonlinear susceptibility X'^' 1 > 10 l6 rrr/V 2 . 

Such predictions are confirmed by detailed numerical solutions of the quantum master equation (5) based on the modeling 
above. In order to take into account the realistic implementation of a quantum nanoplasmonic device, we assume an experimental 
configuration allowing for a coherent driving of the localized plasmon, which can be described by the Hamiltonian 

H = Ho + Fe~ iC0pt p* +F* e iC0pt p, (9) 
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Figure 4. Ratio between the nonlinear coupling energy and the plasmon linewidth, U /(/iyb), as a function of the third-order 
nonlinear susceptibility of the infiltrated nonlinear medium, for some structures of Fig. 3. Gap sizes are: 15 nm (D80), 10 nm 
(B60), and 5 nm (B40). The horizontal dotted line corresponds to the value U/(Jiyo) = \ ■ 


where F(t)/h represents the effective rate of plasmon excitation at the fixed external frequency co p (e.g., imposed by a pumping 
laser). As a figure of merit for quantum nonlinear behavior, we focus on the degree of antibunching in the second-order 
plasmon correlations 48 , G^ 2 \t,t') = (t')p(t')p(t)). In Fig. 5(a) we plot the normalized function at zero time delay 

(t = t' — t = 0), defined as g^(0) = G^ 2 \0) / ((p* p)) 2 , under continuous wave excitation (F/h = 0.01 yo) and as a function of 
the nonlinear susceptibility of the infiltrated nonlinear material, which is supposed to have the same refractive index «b = 1 -5 of 
the background (see Supplementary Information). In all cases, the condition g^(0) < 0.5 can be considered as a threshold for 
single-plasmon blockade, in analogy to photon counting statistics in CQED experiments 1 '■ 12,14 . Ideally, when (0) —> 0 the 
probability of detecting two photons at the output of the device is negligible soon after having detected one. This corresponds to 
the single-plasmon-blockade regime, which we have schematically shown in Fig. 1(b). We assume that the external driving field 
is in resonance with the bare plasmonic mode. Indeed, this corresponds to the condition of maximum plasmon antibunching 8 . 
In general, increasing the detuning between the pump and the plasmonic resonance is detrimental for the observation of the 
blockade effect, especially for the case of positive detuning ( co p > ©o), which could result in the direct excitation of the 
two-plasmon state blue-shifted by the nonlinear interaction. 

In agreement with the more qualitative analysis of Fig. 4, values of the order of X' 3 ' ~ 10~ 16 m 2 /V 2 are predicted from the 
results of Fig. 5(a) to be sufficient to reach the single-plasmon blockade regime in the structure with the tightest electromagnetic 
field confinement, i.e. the bowtie structure with a 5-nm gap between the triangle tips. Such nonlinear values can be achieved, 
e.g., in glasses doped with metallic nanoparticles 49,50 . In such cases, care must be taken in the choice of metallic doping, 
since the plasmon excitations of the nanoparticles embedded in the glass matrix might interact with the target surface plasmon 
excitation of the nanodimer, thus making the whole analysis more complicated. In any case, stronger nonlinearities are necessary 
for structures with inter-particle gaps larger than 10 nm. On the other hand, such high values of the third-order nonlinearity 
might be within reach of available materials. For instance, values on the order of ~ 10 13 — 10 4 3 m 2 /V 2 have been 
measured in organic dye molecules in the visible range 51,52 . Therefore, the use of organic molecules appears a promising 
route for the experimental observation of plasmon blockade. In order to avoid complications due to strong interaction between 
surface plasmons and optical excitations of the nonlinear medium, as well as avoiding nonlocality of the nonlinear response, it 
is important to select a frequency range characterized by a significant third-order nonlinearity together with an essentially flat 
first-order response. 

These results suggest that the single-plasmon blockade regime might be within reach. However, experimentally detecting 
this effect might still be challenging due to the extremely short plasmon lifetime. For a decay rate fryo ~ 100 meV, one expects 
T p = 1 /yo ~ 6.6 fs: this is beyond any possible resolution of single-counting detection. However, such an experiment could 
still be performed under pulsed excitation, where sufficiently short pulses would overcome the issue of short lifetime under 
cw excitation. In Fig. 5(b) we show that this is indeed possible: a train of gaussian pulses with 10 fs duration is sent on the 
structure B60 assuming x 3 ~ 10 15 m 2 /V 2 [see Fig. 5(a)]. The suppression of the zero-time delay peak in the unnormalized 
function G {1) (t = t’ — t), calculated with the master equation after a two-time super-operator evolution 9,9 (see Supplementary 
Information for details), is the signature of single-plasmon blockade in this system (in this simulation, t = 12 fs). In particular, 
this device demonstrates a single-plasmon source on demand, since each pulse triggers the re-emission of a single-plasmon 
from the system before a second one can be excited. Given the short plasmon lifetimes, an ultra-high emission rate can be 
achieved in principle (about 5 THz in this case), beyond state-of-art demonstrations based on plasmon enhancement of single 
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Figure 5. (a) Second-order correlation function at zero time delay as a function of the third-order nonlinear susceptibility 
magnitude, for the same structures as Fig. 4. A schematic representation of the proposed experiment is given in the inset, (b) 
Second-order correlation under pulsed excitation for the B60 structure and % 3! = 10~ 15 m 2 /V 2 (corresponding to the circled 
result in the previous panel). The pulse sequence is represented in the top of the figure. 


emitters spontaneous emission rate 51 . As such, this system can be turned into an ultra-fast single-plasmon source conditioned 
to the availability on a seemingly fast laser source. For instance, ultra-high repetition rate femtosecond lasers have been shown 
up to 200 GHz range 54 . We notice that the usefulness of such an ultra-fast source could be currently hindered by the lack 
of sufficiently fast single-photon detectors, although a ps-time resolution for photon pair correlation measurements has been 
shown by use of a streak camera-based technique in the visible/near-infrared range 55 . On the other hand, such a fast system 
response might be employed, e.g., in engineering novel quantum devices, such as single-photon switches or ultra-fast single 
photon detectors, which we leave for further exploration. 

A possible experimental configuration is schematically described in the inset of Fig. 5(a): a waveguide below the nanostructured 
dimer excites the localized plasmon mode, and scattered radiation is collected from a near-field tip (e.g., a scanning near¬ 
field optical microscope probe), where the output is sent to a Hanbury Brown-Twiss set-up for correlation measurements. 
Coincidences counts from the outgoing pulses would reveal the plasmon blockade regime, similarly to experiments performed 
with CQED systems 12 - 14 . 

In summary, we have theoretically shown the potential interest of metallic nanodimers interacting with ordinary nonlinear 
media to achieve the quantum plasmonic regime, as defined by the true interaction between two plasmonic quanta confined 
in the same spatial region. The main contributions given in this work are twofold: on the one hand, we have provided an 
effective method to estimate the single-plasmonic confinement in the hot-spot of the dimer nanostructures. On the other hand, 
we have shown how the single-plasmon blockade can be evidenced by time-resolved second-order correlation measurements 
under pulsed excitation, thus overcoming the intrinsically short plasmon lifetimes, which we believe might stimulate further- 
experimental research in quantum plasmonics, turning short lifetime into a resource. 
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COMPUTING QUASINORMAL MODES 

We hereby give the detailed procedure employed to compute the characteristics of localized surface plasmons in 
metallic dimers. Surface plasmons are identified with quasi-normal modes (QNMs), i.e., the solutions of a non- 
Hermitian electromagnetic problem, which are calculated in the context of the boundary element method (BEM) 
described in Ref. 1. All calculations are performed within the MNPBEM toolbox [2, 3], which constitutes a freely 
available implementation of the BEM. 

According to the BEM formalism of Ref. 1, space is subdivided in several regions Vj (j = 1,2,...) with locally 
homogeneous dielectric functions (in our case, the inner region of the dimer and the external space). The electric field 
in each region is decomposed into an incident and a scattered part, the latter being identified with the field generated 
by (fictitious) charge and current distributions Oj(r) and h,(r) on the enclosing surface, as follows: 

E sc (r,w)=i— / Gj (r — s)hj(s) ds — V / Gj(r — s)aj(s) ds. (1) 

c JdClj Jdnj 

Gj( r) = exp(ifcj|r|)/|r| is the scalar Green function for region j. Distributions cr,(r) and hj(r) are the unknowns of 
the problem, and they must be computed from boundary conditions at the separating interfaces. 

After the surface is discretized into a set of N simpler polygons (triangles and quadrilaterals), surface charge and 
current distributions become A-dimensional vectors. As illustrated in the text, the QNMs of the system are obtained 
by solving the nonlinear eigenvalue problem [4] 


E(w) x = 0 


( 2 ) 


for the (complex) eigenvalue Cj and the right eigenvector x. The matrix E is defined after Eq. (22) of Ref. 1. Once the 
eigenvector x is known, the corresponding charge and current distributions are obtained from the following expressions: 

oj = GJ 1 x, 

hj = i^-hGj 1 A~ 1 (L 1 - L 2 )i, ^ 

where j = 1,2 indicates the inner and outer regions of the particle, respectively, n is the normal vector to each surface 
element, and we still refer to Ref. 1 for the definition of the matrices Gj, Lj, and A. The quasinormal field £( r) is 
calculated (up to a normalization constant that will be fixed in the following) by replacing the charge and current 
distributions of Eqs. (3) into Eq. (1). 

We have solved the nonlinear eigenproblem (2) with the two-sided Rayleigh functional iterative algorithm [4], which 
allows to compute the eigenvalue and the right and left eigenvectors, indicated as x and y, respectively. Starting 
from an initial triplet (xoi ^o, yo)j at each iteration k the values of the three quantities is updated with respect to the 
previous estimate by solving the equations: 


E(wfc)x fc+ i = E'(w fc )x fc ; 
[E(w fc )]*y fc+ i = [E'(a>fc)]*y fc ; 


Wfc+l = Wfc — 


yjfc+i s(£fc)x fc+ i 
yfc+i E , (w fc )x fc+ i' 


(4) 

(5) 

( 6 ) 


The notation E' indicates the derivative dE/dw, which can be approximated by finite differences by storing the values 
of the previous iteration, whereas E* indicates the conjugate transpose. In order to speed up convergence, the initial 
guess xq can be extracted from the solution of the electromagnetic problem with a suitable external excitation (e.g., 



2 


a dipole emitter) near resonance. In addition, the vectors x and y can be renormalized at each iteration to reduce 
overflow errors. As the eigenfrequency w is typically complex, we had to modify the MNPBEM toolbox to extend the 
refractive index of the particle to the complex domain, by means of a first-order analytical continuation of the form 


Time t(w' + lOj") 


n ex p(w') + ioo 1 


, dn ex p (oj ) 
cko' 


(7) 


where n exp is a cubic spline interpolation of experimental data and the derivative of the resulting piecewise polynomial 
is calculated analytically. As stated in the text, we employed the tabulation of the experimental dielectric function of 
evaporated gold published in Ref. 5. 

In order to ensure the proper normalization of the QNM, we follow an implicit approach involving an additional 
dipole source, similar to that presented in Ref. 6. We add a point dipole source with momentum p = px at the 
position r 0 in the center of the dimer gap (see the inset in Fig. 2 of the main text) and we calculate the field scattered 
back by the particle to the dipole for w —> Co. The normalization condition, reported by Eq. (2) of the main text, 
reads as follows: 


E sc (r, ui) 


—co 


p_£M 

2eq(w — w) 


£(r). 


( 8 ) 


Notice that both terms are proportional to the momentum p. The dependence in the left-hand term is implicit. 
The condition can be simplified by replacing the matrix E _1 in the calculation of the scattered field with the polar 
approximation for ui ~ Co [4] 


s _1 M 


l _ xy* 

io — uj y*[dT,/dio]x' 


(9) 


where x and y are the right and left eigenvectors, respectively. This allows to factor out the singular term from both 
sides of Eq. (8) and to recast the normalization condition in the form 


P ‘ £( r o) 


2e 0 1 

Co y*[dE/dw] 




-£ M e 

+i—n ■ [(Li 
c 


L 2 )A “V + iknL 1( j> e ) + (L 2 A _1 E 1 - L x A“ 1 E 2 )A e ]} . 


( 10 ) 


The term in curly brackets is the same as that in the right-hand side of Eq. (22) in Ref. 1 (where the definition of all 
included quantities can be found), and it is a source term constructed from the held emitted by the dipole. 

Equation (8) does not require any further calculations than those already performed in order to solve the nonlinear 
eigenproblem, greatly improving the computational efficiency. We employed surface meshes with up to about 6000 
elements. Typically, a few iterations of the algorithm are enough to converge to the solution of the eigenproblem. The 
computational time for the calculation of a single QNM (including normalization) is of the order of a few minutes 
with a quad-core 3.1 GHz personal computer. 

The MNPBEM toolbox has been employed also to calculate the semiclassical decay rate of the dipole emitter of 
Fig. 2(a) of the main text. The rate is computed through the formula [7]: 

T(w) = TfreeM + 7-Im [p • E sc (r 0 )], (11) 


where E sc (r 0 ) is the field scattered back by the nanodimer at the dipole position. 


DISCRETE DIPOLE APPROXIMATION 

The modification of the dipole emitter total decay rate induced by the presence of the dimer has also been quan¬ 
titatively assessed in the framework of the discrete dipole approximation (DDA) full-wave simulation method, which 
describes the scatterer as an array of polarizable points organized on a regular cubic grid [8]. This method yields 
solutions for the electromagnetic field in response to an incident electric field in the frequency domain, including 
retardation effects. The polarization of each element internal to the scatterer is seen as the result of the interaction 
with the local electromagnetic field produced by all the other dipoles plus the external field of the emitter. By simply 
hxing the dipole position in space, the changes induced by the dimer antenna on the local field experienced by the 
dipole can be rigorously calculated. The validity of the method been largely tested in recent works [9]. 



3 


The simulations for the present work have been performed by using the parallel code ADDA in its final version, 
which implements dipolar sources [10]. The inter-dipole distance is fixed to 0.25 nm - a value that is small enough to 
achieve numerical convergence - and the gold refractive index is taken from Ref. 5, similarly to the BEM calculations 
described before. 


QUANTIZATION PROCEDURE FOR QUASINORMAL MODES 

In the context of the system-batli approach for the quantization of the electromagnetic field in the presence of 
dispersive and dissipative media, the electric field operator is expanded on a continuum of “true” normal modes 
including the effect of the dispersive media, in the following form [11, 12]: 



where f(r, cu) are (vector) destruction operators satisfying the bosonic commutation relation [/j(r, w),/^(r', w')] = 
<5(r — r ')S(uj — oj')5j k (, j , k = 1,2,3), and “h.c.” indicates the Hermitian conjugate operator. The quantity G(r, r', ui) 
is the dyadic Green function of the (classical) electric held, satisfying V x V x G(r, r',cu) — ure(o;) G(r, r', oj)/c 2 = 
IS(r— r'). Each tensor component Gj k is defined as the j-component of the electric held at point r calculated in the 
presence of the dispersive media for a point-like current source at point r' and of the form: J(r) = — z<5(r— r')~k k /{uipo). 
The source corresponds to an oscillating dipole with momentum p = x k /(uj 2 /j, o); as a consequence, we can use the 
result in Eq. (8) to obtain the expansion of the dyadic Green function in proximity to a QNM: 

2 

Gjk( r,r',w) = - — ^£j(r)€ k (r') + G othel .(r, r', cu). (13) 

Zj\UJ — CO) CO 

The term Gother^, r', u>) includes the free-space contribution and those of the other QNMs of the system. Since, 
for ui ~ Re (a)) the dominant term is the pole corresponding to the QNM, we will assume G 0 ther(r, r', u) ~ 0. This 
assumption is further confirmed by the comparison of the decay rates in Fig. 2 of the main text. 

Along the lines of the approach presented in Appendix A.2 of Ref. 12, we introduce the collective modes p(ui) 
according to the following relation: 


r,w)p(w) = \ — [ dr' i/lm^r', w) G(r, r', w)f (r', w). 

V ttso c . J 


(14) 


The modes form a continuum depending only on frequency and satisfy the bosonic commutation relation [p(w), pl(u/)] = 
S(uj — ui’). Moreover, we assume the vectors g(r, ui) to be real. These assumptions will be justified at the end of the 
procedure. Following the calculations in Ref. 12, by replacing Eq. (14) into the commutation relation, we obtain: 


h oP 1 

9j{r,u)g k ( r',w) = ——jImG ii: (r,r». 

7 T£q C 

Then, by further replacing the expansion of Eq. (13) into the latter expression, we arrive at 


gj(r,u))g k (r\u}) = --^-Im 
27T£o 


'£j(v)£ k (v'Y 

l huiQ 

7o 

CO — Cb 

2tt 2 eq 

_(w - w o) 2 + 7o/ 4 . 


£j{r)£ k {r') 


(15) 


(16) 


[Cb = ujq — *7o/2). As a result, indicating as fl 2 (w) the quantity in square brackets in the last term, we can write 


g(r,w) 


/ fcup g(u>) 

V 2 £q 




(17) 


which, together with Ec^s. (12) and (14), gives Eq. (3) in the main text. In Eq. (16) we suppose that the QNM field 
£(r) can be taken approximately as real, as we verified to be the case for all the systems under consideration in this 
work (Ref Im£ near the gap region). In other systems, this approximation might not hold. This latter situation 
corresponds to the case of a non-Lorentzian density of states discussed in Ref. 13, which leads to a pathological master 
equation not in the Lindblad form considered here. 
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CALCULATING THE SECOND-ORDER CORRELATIONS 


For a given quantized field described by destruction (creation) operators A (A'), the figure of merit quantifying the 
single-photon sensitivity is the time-ordered second-order autocorrelation function, defined as [14] 


G?\t,t')=(A\t)A\t')A{t')A(t)), 


(18) 


where t' — t = r > 0, which is normalized as 


7 ( 2 ) 


(f,0 




(19) 


In the manuscript, we have provided results for the numerical calculations of the latter quantities as follows. The 
master equation in Eq. (5) of the main text is recast in the form 

Hr ot] + \ [2 i/3i f - tiAp - pi f i] , (20) 

where the Hamiltonian is rotated with respect to the driving laser frequency and it is represented as 

H ro t = hA A^A + U (i f ) 2 i 2 + F(f)il + F*(t)A , (21) 


with A = ojq — The Hilbert space is truncated to 11 plasmons, largely sufficient for convergence under the weak 
pump amplitudes assumed in this work (Fg/h = O.OI70). 

Under continuous wave (cw) excitation, F(t) = To, the steady state zero-time delay second order correlation can be 
easily calculated as 

<? (2) (0) = ((it) 2 i 2 )/(iti) 2 = Tr{(i f ) 2 A 2 p ss }/n 2 , (22) 

where n = Tr{yUA/3 ss }, and p ss is the steady state solution corresponding to d/5/di = 0. 

Under pulsed excitation, e.g. for a train of Gaussian pulses described by F(t) = Fbexp{ — (f — nTo) 2 /AT 2 } (with 
n = 0, ±1,...), where To and AT are the pulse separation and width, respectively, Eq. (18) is evaluated numerically 
through a superoperator evolution [15-17] 

G^(t,t') = TA{Au t AMt)ti]ti} , (23) 


where indicates the evolution from t to t' with Eq. 20 when assuming the starting time operator as Ap(t)A^. 
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